Linear Magnetoelectric Effect by Orbital Magnetism 
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We use symmetry analysis and first principles calculations to show that the linear magnetoelectric 
effect can originate from the response of orbital magnetic moments to the polar distortions induced 
by an applied electric field. Using LiFePCU as a model compound we show that spin-orbit coupling 
partially lifts the quenching of the 3d orbitals and causes small orbital magnetic moments ~ 
0.3/jLb) parallel to the spins of the Fe 2+ ions. An applied electric field E modifies the size of these 
orbital magnetic moments inducing a net magnetization linear in E. 



PACS numbers: 75.85.+t, 71.15.Mb, 75.47.Lx 

The last decade has seen increasing interest in the 
study of coupling between electric polarization and in- 
trinsic magnetic moments in materials [1]. Such mag- 
netoelectric coupling manifests in numerous macroscopic 
phenomena: Two well known examples are so-called 
type-II multiferroism [2] in which the onset of mag- 
netic order induces a spontaneous polarization, and lin- 
ear magnetoelectricity, where an applied electric field E 
(magnetic field, H), induces a magnetization Mj = a^Ei 
(polarization, Pi = ctijHj). Although the two phenom- 
ena are non-reciprocal, that is many multiferroics do not 
show a linear magnetoelectric effect and vice- versa, they 
are believed to share closely-related microscopic mecha- 
nisms. 

First-principles computations have been particularly 
informative in resolving quantitatively the various micro- 
scopic contributions to magnetoelectric response [3-5]. 
The first study [3] extracted the "ionic spin" contribution 
to a, by calculating the change in spin canting caused 
by an E-induced polar distortion [6]. Subsequently, the 
methodology to calculate the "electronic spin" compo- 
nent was implemented, through calculating the electric 
polarization induced by an applied Zeeman H field that 
couples only to the spin component of the magnetiza- 
tion [4]. In this method, the electronic spin response is 
obtained by "clamping" the ions during the calculation; 
relaxing the ionic positions in response to the H field 
yields both the ionic and electronic spin components. In- 
terestingly, and perhaps surprisingly, this study showed 
that the ionic and electronic contributions to a can have 
similar magnitudes. 

These spin-based contributions to a have been shown 
to capture much of the experimental response. For the 
case when the magnetic field is applied perpendicular to 
the spins in a collinear antiferromagnet, the magnetoelec- 
tric coupling, a± is relativistic in origin, resulting e.g. 
from the electric-field dependence of the antisymmetric 
Dzyaloshinskii-Moriya exchange [5, 7]. The calculated 
zero kelvin polarizations are consistent with experimen- 
tal values [3], and the temperature evolution of a± fol- 



lows that of the antiferromagnetic order parameter [7]. 
The behavior of an - obtained when the magnetic field 
is applied parallel to the spins - is more complicated. In 
this case, the Heisenberg exchange interactions between 
the spins induce an electric polarization at finite tem- 
perature which is approximately an order of magnitude 
larger than that from the anisotropic exchange interac- 
tions of relativistic origin responsible for aj_[5]. It has 
been shown that responses calculated within this Heisen- 
berg exchange model [5] agree closely with experiment in 
the region close to Tjy (Fig. 1(a)) [8]. One experimentally 
observed feature is lacking, however: While Heisenberg 
exchange predicts ay — » for T — ^ K, consistent with 
the vanishing parallel spin susceptibility at zero kelvin, 
many magnetoelectrics with collinear antiferomagnetism 
have non-zero an at zero kelvin, and instead follow the 
temperature dependence sketched in Fig. 1(a) (solid line). 
An obvious candidate for the discrepancy is the neglect 
of orbital contributions [9]. 

While the neglect of orbital magnetism in the above 
methods is partially justified by the strong quenching 
of 3d orbital moments which usually occurs in transi- 
tion metal oxides, spin-orbit coupling, H so = XL • S, 
can of course reduce the quenching, and allow a non- 
negligible orbital magnetization. This scenario is likely in 
the collinear antiferromagnet LiFeP04 and in LiCoP04. 
Both these compounds have a substantially non-zero an 
as T — )> and an anomalously large anisotropy of the 
magnetic g-tensor [10, 11]. 

Calculation of the orbital contribution to the magne- 
toelectric response is not straightforward, and only a few 
examples, for limited cases and specific approximations, 
exist in the literature. An early study of LiCoPC^ calcu- 
lated the "electronic orbital" (clamped ion) contribution 
analytically, by determining the change in g-factor with 
electric field using perturbation theory within a single- 
ion Hamiltonian [12]. While giving a non-zero value for 
an at T = 0, this method underestimated its magnitude. 
More recently first-principles finite-electric-field methods 
were used to calculate the electronic orbital contributions 
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FIG. 1. a) Qualitative sketches of the temperature depen- 
dence of a || in collinear antiferromagnetic magnetoelectrics 
such as LiFePCU (solid line) or Cr2C>3 (analogous to this curve 
but with negative zero temperature value) and that calculated 
within a spin-exchange striction mechanism (dashed line), (b) 
The orthorhombic unit cell of LiFePCU contains four Fe 2+ 
magnetic cations (brown spheres) which are coordinated by 
distorted oxygen (red spheres) octahedra. Li and P ions are 
represented, respectively, by green and purple spheres. The 
blue numbers label the magnetic sublattices. The arrows in- 
dicate the screw rotation axis parallel to b and c while the 
black dot indicates the center of inversion. 

to the trace of the magnetoelectric tensor - the Chern- 
Simons term - for C^Os and BiFeC>3 [13, 14]. This con- 
tribution was shown to be negligible with respect to the 
spin contribution in both cases. In this letter we explore 
the remaining "ionic orbital" contribution to the magne- 
toelectric response by calculating the dependence of the 
local, on-site orbital magnetic moments on polar lattice 
distortions using density functional theory [15]. Using 
magnetoelectric LiFeP04 as a model compound, we show 
that this ionic orbital contribution to a is unexpectedly 
large and can explain the anomalous low-temperature be- 
havior observed in certain components of a that were 
previously not understood. 

LiFeP04 is orthorhombic (space group Prima) and its 
unit cell (see Fig. 1(b)) contains four magnetic sublat- 
tices occupied by Fe 2+ (S=2) ions. Each magnetic ion is 
surrounded by strongly distorted polar oxygen octahedra 
for which the only remaining local symmetry is a mir- 
ror transformation perpendicular to the crystallographic 
axis b giving local C s symmetry. At temperatures below 
Tat ~ 50 K the Fe 2+ magnetic moments order in the an- 
tiferromagnetic collinear structure with order parameter 
G = 011—012+013—1114 where m 7 is the magnetization of 
the z-th sublattice. The spin orientation in the antiferro- 
magnetic state is still slightly controversial. Early elastic 
neutron scattering and X-ray diffraction data suggested 
that the magnetic moments are fully oriented along the 
b direction [16, 17]. However, recent neutron scattering 
measurements [18] provide evidence for a magnetic struc- 
ture in which G is slightly rotated from b. In this paper, 
we study only those components allowed with G || b; 
G a ^ or G c ^ would give rise to additional non-zero 



components of the magnetoelectric tensor [19] that have 
not yet been reported. 

The onset of the antiferromagnetic order breaks inversion 
symmetry and allows for linear magnetoelectric couplings 
in the free energy 

$H = X\\G b E a H b and $ ± = \ ± G b E b H a , (1) 

where = c^/G 6 and the subscript denotes whether 
the magnetic field is longitudinal or transverse to the 
collinear magnetic moments. 

a II follows the typical form discussed previously and 
sketched in Figure 1 (a): Decreasing the temperature 
from T/v, ay rapidly increases and reaches a maximum 
at T max « 45 K. Below 

Tmax, ol\\ decreases until 20 K 
at which it becomes almost temperature independent 
with a value of ~ 2 ps/m [20]. Importantly, it does 
not approach zero as T — >• OK. a± has the simpler 
temperature dependence mentioned earlier, increasing 
with decreasing temperature below T/v to reach a 
roughly constant value below 25 K (4 ps/m) [21]. 

We focus on the microscopic couplings which can in- 
duce ay. Phenomenologically, exchange-striction cou- 
plings between electric polarization and spins are al- 
lowed by symmetry and give rise to the term: P a oc 
(oil • rei3 — 012 • 014) (see Tab. I). This coupling results 
in a temperature behavior of an similar to that discussed 
above for Cr 2 03 [5]. We note that the local symmetry 
C s of the crystal field around each Fe 2+ ion has only 
one-dimensional irreducible representations and, there- 
fore, the d orbitals are non degenerate. When the or- 
bital moments are fully quenched the magnetic moment 
at the z-th site is proportional to the spin = 2/i^S^. 
As discussed above, at T = the spins in a uniaxial an- 
tiferromagnet are not modified by Hy weaker than the 
magnetic field necessary to flop the spins. Therefore, the 
electric polarization generated at T = by the above 
couplings in response to Hn is zero. 

Next we analyze the orbital contribution to an. We 
begin by discussing the orientation and size of orbital 
moments in zero applied field. From an atomistic per- 
spective, when H so = XL • S is considered the orbital 
moments are partially unquenched and the magnetic mo- 
ment at site i is: 

m!t=n B (2S? + Ut)=t iB grsr (2) 

where L$ and g^ v are, respectively, the orbital momen- 
tum operator and the gyromagnetic tensor at site z, 
/i, v — a, 6, c and summation over repeated indexes is im- 
plied. For an ion with non-degenerate ground state first- 
order corrections in A lead to g^ v — (2 - AAf) where 
Af u = V hM^l^rM^hM , Here ^ is the ground 

t — 'it €ri Co 

state wave function and e n and ip n are, respectively, the 
energy and the wave function of the n-th excited state 
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TABLE I. Transformation of the four magnetic sublattices 
(second to fifth column) under the three generators of the 
space group (modulo a primitive translation) of LiFePO^ 
inversion /, two fold screw rotations around the c axis 2 C , 
and around the b axis 2^ (see Fig. 1). Columns six to eight 
show the transformation of three components of the rank 2 
axial tensor A* at the i-th magnetic sublattice. Here the 
subscripts refer to the change of magnetic sublattice, e.g. 
2 c (A3 b ) = -A^(3) = -Af. The last three columns show 
the transformations of E. 




-0.06-0.04-0.02 0.02 0.04 0.06 
E a (V/A) 
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of the Fe 2+ ion at site i. Since the magnetic moments 
are parallel to b we consider the components A^ b . The 
transformations of these components under the genera- 
tors of the space group (modulo primitive translations) 
are listed in Tab. I, where we see that Af = Af = and 
A bb = A bb at every magnetic sublattice [19]. 

The mean values of the orbital parts of the magnetic 
moments induced by the antiferromagnetic ordering are 
tf L)i = -XAf(S b ). For d 6 ions, A < 0, therefore the 
orbital moment is parallel to the spins in every magnetic 
sublattice and like the spins, gives rise to zero net mag- 
netic moment. 

Next we consider the case E / 0. Electric-field- 
induced polar lattice distortions modify the crystal field 
around each Fe 2+ ion and the energies e n = e n (E). Ex- 
panding Af to first order in E one obtains: Af ,u (E) = 
Mt' v (0)+E a dE«Att''', where 



^0|^|^n>(^n|^1^0)a(Ae, 



(Ae n )2 



dE c 



<T (3) 



and Ae n = e n (E) — eo(E). €p U are the remaining terms 
containing derivatives of wave functions with respect 
to E a . The transformations of the derivatives d E ^A^ b 
under the space group of LiFeP04 can be obtained 
from those of Af 6 and those of E [19] in Tab. I. From 
these transformations we obtain d E aA\ b = d E aA b 3 b = 
-d E aA b 2 b = -dE°Af = d E aA bb . Therefore, the response 
of the average orbital-induced magnetic moment to an 
electric field along a gives rise to a net magnetization 
along b 

M(l) = VBd Ea A bb ((S b ) - (S b ) + (S b ) - (S b ))E a (4) 

that at T = gives fi b (L) = Aji B Sd E a A bb E a . 

To calculate the strength of the linear magnetoelectric 
coupling arising from this mechanism, we perform 
first principles calculations using the Vienna ab initio 
simulation package (VASP) [22]. We use a plane- wave 
basis set for the expansion of the electronic valence wave 



FIG. 2. Calculated electric- field dependence of the net orbital 
magnetic moment per unit cell. E || a (panel a)) results in 
an orbital magnetization along b (ay) while E || b (panel 
b)) produces a net orbital magnetic moment along a (aj_). 
Blue dots and red squares are calculated values at J=l eV 
and J=0 eV respectively, while the straight lines are linear 
fits to the calculated values. The cartoons on the right panels 
show the size and orientation of the orbital magnetic moments 
(gray arrows) of Fe 2+ (brown spheres) when the electric field 
(yellow arrow) is applied along a (panel c) and b (panel d). 
In the cartoon the size of the effect is increased for clarity. 



function and PAW [23] potentials for the treatment of 
core electrons. The exchange-correlation potential is 
described within the local-spin-density approximation 
plus a rotationally invariant Hubbard-^/ (LSDA+£7) 
with a U value of 5 eV, and J values between and 1 
eV. Calculations are performed at the experimental unit 
cell volume of 291 A 3 [17]. We first relax the structure in 
the absence of spin-orbit coupling and then we include 
spin-orbit coupling to calculate the orbital magnetic 
moment. We obtain an orbital moment = 0.306/i# 
parallel to the spins when we use a J value of leV. 
We note that the magnitude of the magnetic moment 
depends on J and on the PAW sphere radius used as 
discussed in Ref. [19]. 

To calculate the ionic orbital response - that is the 
change in orbital magnetic moments when the ions are 
displaced by an applied electric field - we adapt the 
framework introduced in Ref. [3] to obtain the ionic spin 
response. As in Ref. [3], we shift the equilibrium positions 
v t of the ions by Arf = E" o,,, 1 ,.,^,: where 
is the inverse of the force constant matrix after the acous- 



tic modes are traced out and are the Born effective 
charges, both calculated in the absence of spin-orbit cou- 
pling. Since we aim to separate the orbital from the spin 
contribution, we constrain the orientation of the spins 
to lie along the b direction, which we call the "clamped 
spin" approximation (note, however that the magnitude 
of the spin is unconstrained.) After making the Arf 
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distortions from the equilibrium zero-field positions, we 
relax the electronic density with spin-orbit coupling in- 
cluded and calculate the resulting orbital magnetic mo- 
ments. 

Figure 2(a) shows the evolution of the calculated net 
orbital magnetic moment /i^) = J2i=i 4 ^(L)^ of one un it 
cell of LiFeP04 for an electric field applied along a with 
J = 1 eV (blue points) and J = eV (red points). (Note 
that, while the electric field is applied perpendicular to 
the spins, this corresponds to the parallel component of 
a, since the magnetoelectric response is off-diagonal). We 
find that at non-zero electric field the orbital magnetic 
moments remain parallel to the spins, and consistent with 
Eq. (4) the change of their size is opposite for odd and 
even magnetic sublattices giving rise to a net magneti- 
zation. The linear fits of the E a responses of the orbital 
magnetization at J = 1 eV (blue line) and J = eV 
(red line) give ay = 2.3 ps/m and a\\ = 9.3 ps/m respec- 
tively. The an value for J = 1 eV is reasonably close 
to the experimental value of an ~ 2 ps/m at T = 
K [21, 24]. This value of J is consistent with Ref. [25] 
which showed that it is necessary to use J > 0.6 eV to 
obtain the correct magnetic easy axis. To summarize this 
section, we find that the calculated zero kelvin ionic or- 
bital contribution to an has a value which is consistent 
with the measured value of an. We suggest, therefore, 
that the previous discrepancy between the measured zero 
kelvin magnetoelectric response and the calculated spin- 
only response can be explained by this contribution. At 
non-zero temperatures, contributions to ay that are in- 
active in the absence of thermal fluctuations, have to be 
taken into account. These terms comprise the electric 
field dependence of single-ion anisotropy, which has the 
same nature as the orbital magnetic moment, as well as 
the Heisenberg interactions mentioned earlier. 

Finally, we investigate the ionic orbital contribution to 
£*_L> by calculating the effect of an electric field applied 
along b. While the spin-only contribution was not incon- 
sistent with experiment in this case, contributions to a± 
from the electric field dependence of fi^i nave n °t been 
previously investigated and might also play a role. First 
we use similar symmetry arguments as those used for 
an to find constraints on d E bA^ v . From Tab. I we find: 
d Eb Af = d Eb Af = -d Eb Af = -d Eb Af = d Eb A a \ 
d Eb Af = d Eb Af = -d Eb Af = -d Eb Af = d Eb A cb and 
d Eb A bb = 0. On one hand, we note that the transforma- 
tion properties of d Eb Af b are identical to those of d E aA bb . 
This allows for a linear dependence of the orbital magne- 
tization along a when the electric field is applied along 
b: fi a = 4/j jB E b d Eb A ab \(S b }\ where \(S b }\ is the absolute 
value of the average spin component along b. In contrast, 
the transformation properties of d E b A^ 6 , together with 
the spin ordering of LiFeP04 show that the change in or- 
bital moment along c under an applied E b field have op- 
posite sign for sublattices 1, 4 compared with 2,3, yield- 
ing zero net moment in this direction. 



To obtain the size of the ionic orbital contribution to a± 
we perform ab initio calculations using the same method 
discussed for an but with the electric field applied along 
b. As before, we adopt the clamped-spin approach, and 
constrain the spins in the b direction. The resulting cal- 
culated values of net orbital magnetic moment are shown 
in Fig. 2(b) as a function of E b . Here blue and red points 
show results for, respectively, J = 1 eV and J = eV. 
Even when the spins are constrained to be parallel to the 
b axes, the applied Ey induces a canting of the orbital 
magnetic moments from the b direction. In agreement 
with the constraints found for d Eb Af b the resulting cant- 
ing is uniform along the a axis for all magnetic sublattices 
giving rise to a net magnetization linear in E^. Further- 
more, as predicted using the transformations of d Eb A^ b 
for finite we observe a tiny staggered canting of the 
orbital moment along c which gives rise to zero net mag- 
netization. The solid lines in Fig. 2(b) are linear interpo- 
lations of the calculated values and give linear magneto- 
electric responses of 1.9 ps/m and 9.7 ps/m for J = 1 eV 
and J = respectively. To these values, which contain 
only the ionic orbital magnetoelectric effect, one should 
add the spin-only contribution to aj_, which in contrast 
to the case of ay does not vanish at T = 0. These in- 
clude the rotation of easy axis anisotropy, that shares the 
same origin as the canting of orbital magnetic moment, 
as well as Dzyaloshinskii-Moriya interaction. Using the 
approach described in Ref. [4] , which includes these con- 
tributions but not the orbital moment part, we obtain 
for J = 1 eV a value for a± of 2.6 ps/m with sign oppo- 
site to the orbital one. Importantly, these considerations 
can also be used to describe resonant excitation of waves 
of oscillating magnetization M || a with an oscillating 
electric field of a light wave E \\ b : resulting in the so- 
called "electromagnon" peaks in optical absorption [26]. 
Thus the coupling between the orbital magnetic moment 
and electric field gives rise to both static and dynamic 
magnetoelectric effects. 

In summary, we have shown that a linear magneto- 
electric effect can arise from the dependence of orbital 
magnetic moments on the polar distortions induced by 
an applied electric field, the so-called "ionic orbital" con- 
tribution to the magnetoelectric response. We presented 
a symmetry analysis which allows the components of 
a^ v for which this effect exists to be determined, and 
a methodology which can be used to calculate ab initio 
those components at T = 0. We applied the methodol- 
ogy to LiFeP04 and resolved the previous discrepancy 
between previous calculations of the spin-only contribu- 
tions and experiment for ay. Our results show that the 
orbital contributions to the magnetoelectric response can 
be comparable in size to the spin contributions of either 
relativistic or exchange-striction origin in 3d transition 
metal compounds. As suggested by Eq. (4), the temper- 
ature dependence of the magnetoelectric effect caused by 
orbital magnetism coincides with that of the order pa- 
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rameter which, added to the temperature dependence of 
magnetoelectric effect originated by striction gives rise 
to a qualitative agreement for various collinear antiferro- 
magnets such as Cr 2 3 [27], LiCoP0 4 [28] and TbP0 4 
[29]. 

Furthermore, we note that if such coupling between 
orbital magnetization and polar distortion is allowed 
by symmetry, its strength does not depend solely on 
the strength of the spin-orbit interaction. As shown 
in Eq. (3), from a single ion perspective, the strength 
of such an effect is determined by the energy gap be- 
tween the ground state and the excited states for which 
(^ol^lVVi) 7^ and also by the dependence of the en- 
ergies of ionic orbitals on polar distortions of the crys- 
tal field. This suggests that large magnetoelectric effect 
due to the orbital moment correlates with the enhanced 
anisotropic g-tensor and the anisotropy of the magnetic 
susceptibility in the paramagnetic state. In particular, 
large response of orbital magnetism to an applied elec- 
tric field might be found in compounds with reasonably 
small electronic gap, containing magnetic ions with large 
spin-orbit coupling and with low symmetry polar oxygen 
coordination. 
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Effect of small rotations of G away from b. The trans- 
formations of the magnetic field (H), the order param- 
eter (G) and the electric field (E) under the generators 
of Prima (modulo primitive translations) are shown in 
Tab. SI and in Tab. I in the text. The invariant terms in 
the free energy, &me, leading to linear magnetoelectric 
couplings can be written as 

§me = G a (\\H a E a + \\ 1 H b E b + X a ±2 H c E c ) 
+G\\\H b E a + \\H a E b ) 
+G c (\\H c E a + \^H a E c ), (SI) 

where Ay and A^ are functions of temperature. From 
Eq. (SI) it is clear that different components of the mag- 
netoelectric tensor correspond to different components 
of L along the crystallographic axes. Therefore, a small 
rotation of G from the b direction does not affect the 
analysis of ay and a± performed in the text with the as- 
sumption of G || b. 
Symmetry analysis of Af b and d E ^Af b extended to three 
electric field orientations. The effect of the generators of 
Prima on A^ v can be obtained considering the transfor- 
mations of the magnetic sublattices (see Tab. I in the 
text) and the transformations of L. Table SII shows 
the way in which Af changes under these generators 
(modulo primitive translations). These transformations 
imply constraints on some of the A^ ,u components at 
each magnetic sublattice. For instance, the invariance 
under inversion, /, gives: Af = Af and the invariance 
under two- fold rotation around b implies: A J 6 = — A4 6 
which can be both satisfied only for A J 6 = Af = 0. 
Similarly, using Tab. SII, one finds: Af = Af = 0, 

A bb = A bb = A bb = A 66 and A cb = A cb = A cb = A cb = q 

A small applied electric field affects the orbital moments 
by changing A^ v . The expansion of orbital magnetic mo- 
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TABLE SI: Transformation of the magnetic field E and the 
order parameter under the generators of the point group of 
Pnma. 



ment at the i-th sublattice, ^ L ^ v to the linear order in 
E gives 

fi—a,b,c J 

where the second term inside the brackets contains the 
derivatives given in Eq. (3) in the text. 
The relative changes of orbital magnetic moments on dif- 
ferent sublattice due to an applied electric field are de- 

termined by the symmetry properties of d ^ . These are 
listed in Tab. SII and are obtained by combining trans- 
formations of E with those of A f u under the generators of 
Pnma. Similarly to the case of Af 6 , symmetries impose 

constraints on • F° r example, the set of equations: 
d E a Af = -d E aAf and d E * Af = d E *Af, which need to 
be satisfied to ensure, respectively, the invariance under 
inversion and two fold rotation along b (see Tab. SII), im- 
ply that d E aAf = d E aAf = 0. Analogously, we obtain 
that d E aA ab , d E aA cb , d Eb A bb , d E cA ab and d E aA ab have 
to vanish at each magnetic sublattice. Furthermore, the 
non-vanishing terms must satisfy: 

d E aA bb = -d E aA bb = d E aA bb = -d E aAf 

d Eb Af = -d Eb Af = d Eb Af = -d Eb Af 
d Eb Ai b = d E bA 2 b = —d E bA^ b = —d Eb A^ 
d E cA bb = d Ec A bb = -d EC A bb = -d E cA bb . (S3) 

The first two sets of equations in Eq. (S3) together with 
Eq. (S2) and the sign of the components of the spins along 
b at every magnetic sublattice shows that an applied elec- 
tric field along a (b) induces a net orbital moment along 
b(a). 

In the same way, from the last equation in Eq. (S3) one 
can derive that an applied electric field along c affects the 
orbital magnetic moments along b of sublattice 1 and 4 
in the opposite way to those of sublattices 2 and 3 giving 
rise to a zero net magnetization. However, the staggered 
orbital magnetization C (L) = M(l)i-M(l)2~M(l)3+M(l)4 
along b should by symmetry depend linearly on the ap- 
plied electric field. To check this we perform the same 
calculations described in the text but for E || c. The 
obtained electric field dependence of C b for J = 1 eV 
is plotted in Fig. SI (a) and shows that such response is 
approximately double that of the orbital magnetization 
described in the text. In the same way we expect a linear 
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TABLE SII: Transformation of tensor Af b and its derivatives 
with respect to the electric field under the generators of Prima 
(modulo a primitive translation). The superscript indices la- 
bel the coordinate while the subscript index labels the mag- 
netic sublattice. 
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FIG. SI: Electric field dependence of C b L) (a) for J = 1 eV 
(blue dots), (b) J dependence of the size of orbital moments 
at each magnetic sublattice for U = 4 eV (blue dots) and 
U = 5 eV (red dots), (c) Size of the orbital moment as a 
function of integration sphere volume for E — V/A (blue 
dots) and E = 0.06 V/A( green triangles and red squares). 
The lines are guides for the eyes. 



dependence of on E || b. However, such dependence 
is so weak that, as mentioned in the text, the changes in 
orbital moments are above numerical resolution only at 
E = 0.06 V/A. 

We also note that the size of orbital moments strongly 
depends on the values of the on-site Coulomb interac- 
tion U and on the effective on-site exchange parameter 
J. Figure Sl(b) shows the size of calculated orbital mag- 
netic moment at E = as a function of J for U = 4 eV 
(blue dots) and U = 5 eV (red dots). 
Convergence of the orbital moment. Finally, we analyze 
the convergence of the local orbital moments with re- 
spect to the radius of the spheres in which they are eval- 
uated. The modulation of this radius in VASP would 
require the generation of new pseudo-potentials, which 
is not an easy task. Hence we decided to use an all- 
electron method based on linear augmented plane waves 
(LAPW), namely ELK [1], which allows to alter the ra- 
dius of the projection spheres. First, we compared orbital 
moments using both methods and the same radius of the 
spheres finding that in the ground-state structure their 
values agree up to 10 -4 /i#- Next, we varied the radius 
of the spheres in the ELK code. For this test, we used 
a GGA exchange correlation potential to avoid any ad- 
ditional influence of the Hubbard U applied within the 
sphere. Figure Sl(c) shows the obtained value of orbital 
moments as a function of sphere volume at two Fe sub- 
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lattice sites for the case of E = (blue squares) and 
E b = 0.06V/ A (red squares and green triangles). While 
the orbital moment increases approximately linearly with 
sphere radius over the investigated range, the difference 
of the orbital moment is constant, so we can conclude 
that the response is well converged even at the smallest 
sphere radii. It is worth mentioning that the orbital mo- 
ment response obtained using GGA is smaller than that 
from with LDA+£7. This finding may be attributed to a 
volume effect since the structure used was not optimized 
for GGA. 

Additional technical details. All VASP density func- 
tional calculations were performed with an energy cut- 
off E cut = 600 eV. The ionic relaxation, the calculation 
of the force constant matrix and the calculations of the 
orbital magnetism were done using a 2 x 4 x 4 Monkhorst- 
Pack k-point mesh. The Born effective charges were ob- 



tained using both the finite electric field method (after re- 
laxation using 4x4x4 gamma centered k-point mesh) and 
density functional perturbation theory (with a 2 x 4 x 4 
Monkhorst-Pack k-point mesh). The results obtained 
with the two procedures had minimal differences. The 
force constant matrix was obtained using finite displace- 
ments method obtained through the software Phonopy 
[2]. In the LAPW ELK calculation we used a 3 x 5 x 5 
k-point mesh and iterated the density until the change in 
potential was smaller than 10 -6 eV. 
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